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ABSTRACT 

Cosmological simulations of galaxy formation often rely on prescriptions for 
star formation and feedback that depend on halo properties such as halo mass, 
central over-density, and virial temperature. In this paper we address the conver- 
gence of individual halo properties, based on their number of particles N, focusing 
in particular on the mass of halos near the resolution limit of a simulation. While 
it has been established that the halo mass function is sampled on average down 
to N ~ 20 — 30 particles, we show that individual halo properties exhibit signifi- 
cant scatter, and some systematic biases, as one approaches the resolution limit. 
We carry out a series of cosmological simulations using the Gadget2 and Enzo 
codes with N p = 64 3 to N p = 1024 3 total particles, keeping the same large-scale 
structure in the simulation box. We consider boxes of small (Ibox = 8 Mpc h _1 ), 
medium (l box = 64 Mpc h -1 ) and large (lb ox = 512 Mpc h -1 ) size to probe dif- 
ferent halo masses and formation redshifts. We cross-identify dark matter halos 
in boxes at different resolutions and measure the scatter in their properties. The 
uncertainty in the mass of single halos depends on the number of particles (scal- 
ing approximately as iV -1//3 ), but the rarer the density peak, the more robust 
its identification. The virial radius of halos is very stable and can be measured 
without bias for halos with N > 30. In contrast, the average density within a 
sphere containing 25% of the total halo mass is severely underestimated (by more 
than a factor 2) and the halo spin is moderately overestimated for iV < 100. If 
sub-grid physics is implemented upon a cosmological simulation, we recommend 
that rare halos (~ 3a peaks) be resolved with N > 100 particles and common 
halos (~ la peaks) with iV > 400 particles to avoid excessive numerical noise 
and possible systematic biases in the results. 
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Introduction 



Cosmological simulations are a fundamental tool for investigating the fo rmation and evo- 
lution of dark matter halos and of the resulting galaxies ( lBertschingerlll998l ) . With advance- 
ments in both memory and computing power capabilities, simulations can now be performed 
with up to several billion particles. Efforts have been f ocused both on invest i gating the for- 



mation of single halos wi th ultra-high resolution (e.g., iDiemand et al.l 120071 ; ISpringel et al. 



2008at IStadel et al.ll2009l ) and on simulating structure formation in larg e boxes, on the order 



of Gpc 3 , that repr esent the large-scale structure in the Universe, (e.g., ISpringel et al.l 12005 



Yepes et al.l 120081 ) . However, the dynamic range of scales to be resolved is so large, despite 
the continuous hardware improvements, that numerical simulations need to be pushed to 
their limits for many interesting problems. For example, to foll ow the formation of the 
first stars in the Universe, a sub-solar mass resolution is needed (jO'Shea fc Normanl 120071 ; 
Turk et al.ll2009l ). This limits the size of the simulation box to below ~ 1 Mpc 3 , and thus 
large-scale structure effects are missed. At other range of the spectrum, simulations of the 
formation of bright (and thus rare ) high-redshift q uasars observed in all-sky surveys such 
as the Sloan Digital Sky Survey (IFan et al.l 120061 ) require computational volumes > 10 9 
Mpc 3 and resolution sufficie nt to identify the black hole seeds from m etal-free stars, for a 
total dynamic range > 10 12 (ITrenti fc Stiavellil 120071 ; iTrenti et al.ll2008l ). Under such condi- 
tions, information from barely resolved halos needs to be used. In addition, sub-grid physics 
reci pes are often employed , eith e r in the form of semi-analy t ical (post-processing) mo del- 
ing ( jKauffmann fc Charlotl Il998l ; ISomerville &; Primackl Il999r |Pe Lucia &; Blaizotl 120071) or 



as star formation and feedback r ecipes implemented during the run (ISpringel fc Hernquist 
20031 : bppenheimer fc Pavel l2008h . 



Therefore, convergence and validation of the numerical methods used are of fundamen- 
tal importance to establish the reliability of the conclusions drawn from numerical exper- 
iments. Past investigations have addressed two fundamental issues regarding dark matter 
(DM) halos: the convergence of the inner slope of the density profile and the accuracy of 
the halo mass function. The inner-slope problem, aris ing from the absenc e of observational 



evidence of the density cusps predicted in simulations (|Navarro et a 



cus of several studies using extremely high resolution ( IGhigna et al 



2000 



1997 ), has been the fo- 



Power et al. 2003 



Fukushige et al.l 12004 ). but the issu e appears to be sett l ed, with the inner regions expected 
to follow a Einasto density profile ( INavarro et al.l l2010t IStadel et al.l 120091 ). Regarding the 
DM halo mass function, it is now well established that a cosmological simu lation reproduces 
with fid elity the mass function down to halos with N > 20 — 30 partic les (IReed et al.ll2003l ; 



see also Heitmann et al. 



be box-size effects ( jBagla fc Ray 



2006 



Warren et al. 



2005 



2006 



Reed et al. 



excellent agreement with analytical predictions (jSheth fc Tormenl 1 19991 ; I Jenkins et al.ll200ll ). 



2007. 



Lukic et al.l 120071 ) . although there can 



2009). These results are also in 
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As we will demonstrate, the mass function is only an average property of the halo mass 
distribution. A much deeper question pertains to the reliability of properties of single halos, 
derived from a simulation with given resolution, with a goal similar to the study of subhalos 
properties in the Aquarius run ( Springel et al. 2008a ; Springel et al. 2008bl). Individual hal o 
properties affect the variance of the results derived from the simulations (IWarren et al.ll2006l ). 
For example, if semi-analytical formulae and/or star formation recipes are implemented 
in a run, the uncertainty on individual halos propagates to th e derivation of quantities 
such as the star formation rate, the fundamental plane thickness ( iDjorgovski fc Davislll987l ) 
or the tig htness in the relation bet ween the central black hole mass and bulge velocity 
dispersion (IFerrarese fc Merrittl 120001 ) . 



The goal of this paper is to quantify the numerical scatter and identify possible biases 
in the mass of individual halos as a function of the number of particles in the halo. Past 
investigations have characterized the behavior of halo fin ders for idealized syste ms with small 
N, for example by generating discrete realizations o f a iNavarro et al.l ( 19971) profile or b y 



downsampling the resolution of a simulation snapshot ( IWarren et al. 



2006 



iLukic et ajj|2009j ). 



We extend these studies by studying the convergence of halo properties in a fully cosmological 
context, where the convergence properties of the N-body integration are also investigated. 
We carry out a suite of simulations, where higher resolution boxes are constrained to have the 
same phases as low-resolution versions. We cross-identify halos between the different runs 
and measure convergence of their properties, which turns out to be di fferent with respec t 
to the more idealized numerical experiments of I Warren et al.l (120061 ) and lLukic et al.l ( 120091 ) . 
Our study should help the community of numerical cosmologists to quantify the limit at 
which they should trust their simulations, depending on the desired accuracy goal. 

This paper is organized as follows. In Section. [2] we introduce the setup of our suites of 
simulations, whose results are presented in Section. [3J We conclude in Section. H]by discussing 
our recommendations for the minimum resolution required for implementing extra physics 
on dark matter halos. 



2. Numerical setup 



We generate initial conditions using the G raficl package (IBertschingerl l200ll ). with a 
custom modification that allows u s to apply the iHoffman fc Ribaki (119911 ) method over the 
full simulation box and to use the Eisenstein fc Hu ( 19991 ) power-spectrum fitting formula. 
With our customization, we are able to start from a low-resolution version of the initial 
conditions and then refine it to hi gher resolution while k eeping the same large-scale structure. 
We use a WMAP-5 cosmology fcomatsu etaD 12009^ with fi A = 0.72, Vt m = 0.28, tt b = 
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0.0462, erg = 0.817, n s = 0.96 and h = 0.7. We consider three different box sizes: a small 
box of edge lf, ox = 8 Mpc h~ x , a medium box (edge lj, ox = 64 Mpc Yr x ) and a large box 
{hox — 512 Mpc h^ 1 ). The simulations have a range of total particles from N p = 64 3 
to N p = 1024 3 (the number of particles in a single halo is instead indicated as N). The 
small-box simulations start at redshift z = 199, while the medium and large-box simulations 
begin at z = 100. The boxes have periodic boundary conditions. Details on the specific 
simulations, including their mass resolution, are shown in Table [TJ 

As our code of choice we use t he particle-mesh tree code Gadget2 in its "lean" version 
( ISpringelll2005r ISpringel et al.ll2005l ) to carry out dark-matter only simulations. The softening 
parameter is set to e = l^ ox j (42.5N p ), allowing us to achieve a good spatial resolution of virial- 
ized halos with a small number of particles. For comparison, we al so carry out a subs et of the 



runs (the medium-box series) using the hydrodynamic code Enzo ( Brvan et al. 1995T). Enzjf 



uses t he block-structured adaptive mesh refinement (AMR) scheme of iBerger fc Colella 
( 119891 ) to achieve high spatial and temporal resolution, and it combines an N-body adaptive 
particle-mesh solver for dark matter dynamics with a Piecewise Para bolic Method (PPM) 
hydro solver that has b een optimized for cosmo logical applications (jColella Woodward 
1984] ; iBryan et al.lll995l ). Because our primary goal is to test the DM halo dynamics of 
Enzo, we do not include gas in the Enzo simulations. The number of top-grid cells in our 
Enzo runs is equal to the total number of particles. 

The main difference between the Gadget2 and Enzo under these conditions is their 

1 /3 

spatial resolution. The force resolution in Enzo is twice the grid size, that is 2lb ox /N p . In 
Gadget2, the force becomes unsoftened at a distance of about 3e, or at ~ hox/ (14Ay ) in our 
simulations. The Gadget2 runs thus have about 30 times better spatial resolution than their 
Enzo counterparts, if a uniform grid and no AMR is used in Enzo. To investigate the effect 
of force resolution on the properties of halos, we have carried out a subset of the Enzo runs 
(N p = 64 3 to N p = 256 3 ) allowing up to six levels of AMR. This improves the force resolution 
by a factor up to 2 6 , thereby reaching a maximum force accuracy comparable to that attained 
by Gadget2. Enzo still has a lo wer force resolution in regions with overdensities below the 
critical threshold for AMR (see lO'Shea et al.ll2005l ). While many applications of Enzo 



rely on a ggressive use of AMR, for example in the c ontext of the formation of Population 



III stars ( lO'Shea &: Norman! 120071 : iTurk et al.l 12009), a growing number of investigations 



consider runs with uniform resolution (IRegan et al. 



2009 



20071 : IPaschos et al.ll2009l : iTytler et al. 



Norman et al.ll2009l ). Especially for studies of the Lya forest, it has been shown by 



Regan et al.l (120071 ) that disabling AMR provides an order-of-magnitude speedup, while only 



2 http://lca. ucsd.edu/projects/enzo 
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introducing ~ 5% errors in the Lya flux power spectrum. 

We save snapshots of the simulations at regular redshift intervals (z = 6,4 , 1, 0) and we 
identify DM halos with a friends-of- friends (FoF) algorithm ( Davis et al. 1985h with linking 



lengt h 0.2. We also analyze a subset of runs with the Amiga halo finder ( iKnollmann fc Knebe 



20091 ) that includes a boundness check for the halo particles. In addition, the Amiga halo 
finder provides detailed information on each halo, including their density profile and spin 
parameter. 

We consider halos with at least 32 particles. To cross-identify the same halo in two simu- 
lations at different resolution, we match individual particle identification numbers (IDs) that 
ar e representative of the initial particle positions. Our method is similar to that discussed 



in 



Springel et al.l ( I2008af ). The ID of every particle in the lower resolution realization (with 
N = Niow) is used to calculate its corresponding n c = N high /Ni ow "chikT'B particle IDs in the 
high-resolution simulation (with N = N^igh)- From the list of particle IDs in each halo we 
can thus ascertain whether that halo has one or more counterparts in the higher-resolution 
run. Similarly, given a halo in the high-resolution simulation, we can determine the presence 
of any low-resolution counterparts. Note that the relation between halos in snapshots at 
two different resolutions is not necessarily one-to-one or one-to-zero. In fact, multiple halos 
can be the counterparts of a single larger halo, especially in the process of merging. For 
any given low-resolution halo, we identify its high-resolution counterpart by considering the 
high-resolution halo that has the largest number of individual "children" particle matches. 
In passing we note that if matching of halos is based instead on their positions, there is no 
guarantee of either positive or unique identification. This affects especially common low- 
mass objects in the proximity of larger halos, because a change in resolution can lead to 
different tidal forces and changes in the rate of mergers. 



3. Results 

In Figured] we show the scatter of individual halo masses in our medium-box simulation 
at z = when the resolution progressively increases from N p = 64 3 to N p = 512 3 . We plot 
the ratio of low-to-high resolution mass of halos Mi om jM^\2 as a function of the halo mass in 
the highest-resolution run M 512 . As the resolution is increased, it is clear that the halo mass 
is measured with progressively higher accuracy. In Figure [2] we show the scatter as a function 



3 If the N p is increased by 2 3 , then each low-resolution particle has 2 3 counterparts, or "children" at high- 
resolution. The average position and velocity of the "childre n" corresponds app roximately, but not exactly, 
to the position and velocity of their low-resolution "parent" (|Bertschingerll200lh . 
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of the number of particles in the low resolution run (Ni ow ). The median mass of halos is 
correct down to about 100 particles, while it tends to be underestimated for the smallest 
halos. However, the halo mass function remains consistent as shown in Figure [3] (see also 
Reed et al.l 120071 ). We quantify the dimensionless scatter, C(^iow)i around the median for 
halos with Ni ow particles by considering one half of the symmetric la interval that encloses 
from 16% to 84% of the M low /M high points for halos with N low - AN <N< N low + AN. We 
chose AiV > 10 adaptively to ensure that the distribution is well sampled. The scatter around 
the median grows steadily as the number of particle decreases. For halos with less than 10 3 
particles, their mass has > 20% uncertainty at 68% confidence level. The convergence of the 
mean h alo mass down to a small number of particles is consistent with previous resolution 
studies (IReed et al.l 120031 ) that demonstrated that the halo mass function of a simulation is 
correctly sampled down to such low particle number. 

Although the mean halo properties we find are consistent with those reported in earlier 
investigations, the individual scatter of halo masses measured from cos mological simu l ations 
diff ers from the estima tes based on i dealized exper i ments carr ied out by I Warren et al.l (120061 ) 
and lLukic et~aD koO^ \. First, both lWarren et al.l tood ) and lLukic et all ( 120091 ) report that 
FoF halo finders tend to overestimate the mass of discrete realizations of a halo at low N. In 
actual simulations, the opposite behavior is observed. The mass of small N halos is in fact 
underestimated at low resolution. By comparing the N p = 256 3 realization of the medium 
box against the N p = 512 3 realization, we find that the average of Mi ow /M high for halos with 



32 < N < 64 is 0.84 and the median is 0.77. From Table 1 and Equation 3 of IWarren et al. 
(120061 ) we would have instead expected (Mi ow /M high ) > 1.1. We obtain a different result 



because the mass of the halo in a simulati on depends no t only on th e convergenc e prop erties 
of the halo-finder algorithm, explored by I Warren et al.l (120061 ) and iLukic et al.l (120091 ) , but 
also on those of the code that resolves the non-linear gravitational dynamics leading to the 
formation of the halos. A lower number of particles, and thus a reduced force and spatial 
resolution, suppresses high frequency modes of the effective power spectrum of the simulation, 
pro ducing; an underestim ate of the mass of halos near the resolution limit of the run (see 
also lO'Shea et al.l (120051 ) for a similar finding in the context of the comparison between Enzo 
and Gadget2). A second differe nce between our find ings and those publishe d earlier is in the 
amplitude of the scatter. Both I Warren et al.l (120061 ) and ILukic et al.l (120091 ) observe in their 
halo experiments a scatter in the measure of halo ma sses that is a f a ctor 2 lower than what 
is realized in a cosmological simulation. For example, I Warren et al.l ( 120061 ) measure an 18% 
relative l-a error for their synthetic N = 50 halos. In our simulations, the 1-cr relative error 
is > 55% for halos with N ~ 50 as shown in the bottom left panel of Figure [21 Again, this 
is not surprising, because our results are affected primarily by the different resolution in the 
cosmological simulations, rather than by the convergence properties of the halo finder. 



-7- 



Interestingly, the considerable scatter in the individual halo mass is left essentially 
unchanged if the DM halo catalogs are pruned of unbound particles when iV > 100. This 
is shown in the right panels of Figure [2J which are the equivalent of the left panels but 
obtained using halo catalogs from the Amiga halo finder that includes a boundness check for 



membership of particles to a halo ( iKnollmann fc Knebd 120091 ) . Removing unbound particles 



helps only at the lowest end of the resolution for halos with iV < 100; for example ^Amiga(N = 
40) ~ 0.4 while £,fof{N = 40) ~ 0.55. As expected, the overall number of halos identified 
in a snapshot above a given halo mass is slightly lower (~ 5%) when unbound particles 
are removed. For example, in the medium box with N p = 512 3 particles, there are 98,011 
halos with A^ > 32 particles identified by the friends-of- friends halo finder and 93,364 by the 
Amiga halo finder at z = 0. 

The scatter in halo masses remains largely unchanged when we consider different box 
sizes and redshifts, as shown in Figure HI There is a moderate tendency for rare halos to be 
better resolved at a given number of particles compared to their more common counterparts, 
especially when A" < 10 2 (see Figure [5]). This is highlighted by quantifying the rarity of halos 
using the extended Press-Schechter variable v = 5^/a 2 (M). For example, the very common 
(low v) halos in Figure [5] have considerable more scatter at A^ ~ 100 than rarer halos (high 
v) with a similar number of particles. In addition, common halos with ~ 10 2 particles may 
be in reality part of a larger halo when the numerical resolution is increased (see the points 
at Mi ow /Mhi g h ~ in the bottom left panel of Figure [2]). This effect does not happen for 
rarer halos (see the upper left panel of Figure [2]). This is not surprising, because if a halo 
originates from a rare peak, then it is more likely to be the dominant gravitational source in 
its surroundings and the dynamics of its own particles is primarily governed by self-gravity. 
In contrast, more common halos are likely to be surrounded by at least comparably massive 
neighbors, and they might be more affected by tidal-field errors. 

The scatter of individual halo masses is reduced as the number of particles in a halo 
increases, but it remains considerable, even when a halo has A" > 10 3 particles (see Figure [6]). 
The scaling of convergence with the number of particles can be understood with a simple 
analytical model. For the purpose of computing the total mass of a halo, the particles 
more likely to be affected by errors in their dynamics are those initially located at the 
periphery of the halo over-density. In the linear regime at z ^> 1, when the density field is 
quasi-homogeneous, a spherical region that contains A" particles has Nb OU ndary edge particles, 
where: 

f3N\ 2/3 

-^boundary ~ 47T I — J (1) 

V 47r / 

Assuming that a fraction of the edge particles is affected by numerical resolution, then the 
dispersion of the mass of an individual halo scales as £(N) oc N boundary /N ~ A^ 1 / 3 . If 
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we assume that, on average, about half of the edge particles are susceptible to change of 
membership when the resolution is increased, then we expect a ~ 20% uncertainty in the 
mass of a halo with 1000 < N < 2000 particles, in reasonable agreement with the scatter we 
measure (~ 15% in Figure [5]). This means that to reduce the typical uncertainty on a halo 
mass below 10%, iV ~ 5000 particles are required. Because this scaling depends only on the 
surface-to-volume ratio, the iV -1 / 3 trend is predicted even if the spherical assumption for 
the collapse is relaxed and more reali stic models for the formation of halos are considered, 
such as the ellipsoidal collapse model ( iSheth et al.ll200ll ). Figure E] shows that the empirical 
measurements for £ (N) do indeed show that a good fit of the overall distribution is given by: 



( N 
0.15 ( 



1000 



-1/3 



(2) 



From Figures [T] and |2] it can be seen that some halos with a large number of particles 
(N > 10 4 ) can occasionally have a large variation in their mass when the same box is 
resimulated at higher resolution. These are halos in the process of merging, as shown in 
Figure [7] for a halo with 8192 particles in the £& OT = 64Mpc h~ l box (jV p = 256 3 ) at z = 0. 
While the FoF halo finder flags the halo as a single entity at low resolution (there is a bridge 
of particles connecting the two main components), at high resolution the merging is slightly 
delayed, so that the two sub-components are still separate halos. Of course, the opposite 
condition may also be realized, with two individual halos identified at low resolution and a 
single halo at high resolution. Such ambiguity in defining a halo cannot be avoided unless an 
additional diagnostic is used in addition to halo-finding algorithms (such as halo profilers or 
indicators for an irregular morphology). Nevertheless, the fraction of halos in the process of 
undergoing a merger in any single snapshot is small, and the measure of £ (N) is not affected 
by outliers in the distribution of Mi ow /Mhi g h- 

The dimensionless scatter £ we measure appears larger by about a factor two compared 



to the scatter quantified for subhalos in the Aquari us simulation (Figure 1 6 in ISpringel et al. 



2008al : see also the supplementary infor mation in ISpringel et al.l l2008bl ). The better con- 
vergence of subhalo properties found by Springel et ali ( 2008a ) is not surprising. In fact, 
subhalos are the remnants of initially more massive halos that have been stripped of their 
less bound particles, both during the merging with the main halo and by tidal forces, once 
the subhalo is orbiting inside the parent halo. Loosely bound particles are more likely to 
be added or removed from a halo as a result of a change in resolution. In addition, once a 
subhalo is part of a larger halo, its orbit becomes defined by the parent halo potential, and 
further merging with other subhalos is highly unlikely. Scatter in the individual properties of 
subhalos is thus not affected by the ambiguity in defining a halo undergoing a major merger, 
contrary to what might happen for the halos studied here (Figure [7]). 
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The convergence of individual halo properties appears worse in Enzo runs when N < 400. 
Figure [3] shows t h at the halo mass function of such runs deviates from both the analytical 
Sheth fc Tormenl (Il999l ) prediction and the higher-resolution realization when N < 400. 



Increasing the force resolution by switching on AMR helps with respect to Enzo Unigrid 
runs, but there is still a significant number of halos missing at N < 400. This trend identified 
in the halo mass function is also clearly visible at the level of individual halo masses (see 
Figure |8]). The median of the distribution of Mi ow /M high approaches 1 only at iV ~ 1000, in 
sharp contrast to the better convergence properties foun d in the Gadget2 runs. This result 



is no t surprising and has been already previously noted ( jO'Shea et al.ll2005l ; lHallman et al. 



20071 ). The gravity solver in the Enzo code does i n fact suppress sma ll-scale power at very 



high redshift, before AMR refinement is triggered ( jO'Shea et al.ll2005l ). The convergence of 
the halo mass can b e improved by increasing the dimension of the top-level grid, as shown by 



O'Shea et al.l (120051 ) . However, this might not always be possible for the largest runs, when 



the limiting factor of computational resources is often the availability of memory rather than 
processor speed. 

Figure [9] show the convergence of other individual halo properties, obtained by profiling 
the halos with the Amiga halo finder. The core density (P25), defined as the average density 
within the radius containing 25% of the total halo mass, is reported in Figure M for one 
snapshot at z = from the medium box. As expected, the errors on p 2 5 are larger than 
those on the total halo mass, and there is clearly a systematic bias: low-resolution halos with 
small N have a significantly lower p 25 compared to their realization at higher resolution. The 
median of p 2 5 low/ P25 high is as low as ~ 0.3 for iV ~ 32. This bias has a major i mpact if 



the ra te of star formation is scaled from the core density, for example by using the [Schmidt 



( 119591 ) law. The systematic underestimation of core density in low-resolution runs is likely 
due to two-body relaxation at the center of poorly resolved halos. The scatter of the virial 
radius, r vir , is shown in Figure IHJ This quantity appears to be well defined, and has no bias 
and a small scatter down to N = 32. The stability of r„; r is not surprising as r vir cx 
Hence, an error on the mass does not severely affect the associated radius. Figure [9] shows 
the spin of individual halos, measured by the dimensionless parameter A = J\E\ 1 ^ 2 /GM 5 ^ 2 , 
where J is the angular momentum of the halo, E its total energy (gravitational and kinetic) 
and G the gravitational constant. At small N, Xi ow appears systematically higher than Xhigh 
and has a large scatter. This is not surprising: small N systems always have some residual 
angular momentum, even if the particle positions and velocities were to be drawn from the 
distribution function of a non-rotating system. The bias in A for small N may be important 
in semi-analytic modeling of galaxy formation, in case the angular momentum of a halo is 
used to determine the presence of a disk. 
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Conclusions and Discussion 



In this paper we quantify the convergence of individual halo properties, as the resolution 
of a cosmological simulation is increased while maintaining the same large-scale structure of 
the coarser run. We confirm past investigations of the convergence of the global mass function 
down to our resolution limit of N = 32 particles in runs carried out with the particle-mesh 
tree code Gadget2. At the same time, we demonstrate that the scatter in individual halo 
masses measured from our set of cosmological simulations is qualitatively and quantitatively 
different from that reported based on the analysis of the Friends-of Friends halo finder 
on mock halos and on downsampled snapshots of an individual cosmological simulation 
(IWarren et al.ll2006l ; iLukic et al.l 120091 ) . Those experiments highlighted the tendency of the 
FoF halo finder to overestimate the mass of poorly resolved halos. We show instead that, 
in an actual resolution study of cosmological simulations, the mass of low-N halos tends to 
be underestimated. This behavior in the convergence of low N halo properties is determined 
primarily by the finite accuracy o f the gravity i n tegra ti on, rather t 



of the FoF algorithm explored in IWarren et al.l ( 120061); ILukic et al 



r an by the properties 
d2009|) . Jn fact, the 



results obtained in Figure [2] with the Amiga halo finder ( iKnollmann fc Knebdl2009l ) are fully 
consistent with those obtained with the FoF finder, despite the fact that Amiga identifies 
halos from the topology of isodensity contours. 

In addition to the different direction of the systematic bias, the scatter in the individual 
halo masses is more than a factor two higher than reported previously for idealized halos. 
The mass of halos resolved with N ~ 1000 particles shows a ~ 15% scatter when the halos 
are resimulated at higher (8x) resolution. Halos with a smaller number of particles have a 
larger uncertainty in their masses, typically of the order of ~ 30% for iV ~ 100. The la 
relative uncertainty in the mass of a single halo with N particles scales approximately as 
£ ~ iV -1 / 3 (see Eq. [2] and Figure [6]). As the resolution limit of the simulation is approached, 
halo masses show less scatter if the small- iV halos are rare, and thus the simulation box is 
mostly composed of particles still near the regime of linear evolution. The mass of rarer halos 
can instead be measured relatively well, even with a low number of particles (see Fig. [5]). The 
core density and the spin of small- N halos exhibit a systematic bias in addition to significant 
scatter. The virial radius of a halo appears to be the quantity with the smallest scatter, no 
bias down to N = 32. 

We carried out our analysis using a code — Gadget2 — that can efficiently reach a high 
force resolution (a small fraction of the initial inter-particle distance) at low computational 
cost. Simulations done with the AMR code Enzo are prone to a systematic bias in individual 
halo mass when N < 1000. Resorting to Adaptive Mesh Refinement lessens but does not 
resolve this problem (see Figure |H]). Increasing the top-grid dimension helps to improve 
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convergence of individual halo masses with N ~ 100 (see lO'Shea et al.ll2005l ). 



Our investigation highlights the importance of carrying out careful resolution studies 
to validate the conclusions of numerical simulations. This is particularly important when 
sub-grid physics recipes are implemented within a numerical simulation or when analytical 
models of galaxy formation are constructed from the simulation snapshots. If these recipes 
prescribe to populate halos with a small number of particles, one might obtain the correct 
average behavior of the sample (because the halo-mass function converges down to halos 
with ~ 30 particles). However, the limited numerical resolution is likely to introduce extra 
scatter in the properties of the end-products of such simulations. This extra scatter may 
also propagate from low-iV progenitors to descendant halos (and galaxies) of larger mass, for 
instance if the star-formation recipes depend critically on the initial metal enrichment. 

Furthermore, systematic biases are possible. One obvious example is the rate of star 
formation calculated from the central density of a halo. Figure [9] clearly shows that the 
core density of halos with iV < 400 is underestimated by more than 20% (and by more 
than 75% at iV ~ 32). Even recipes based on halo mass alone can lead to systematic 
biases. One example is the ratio of Ly-a to stellar luminosity in a simple model where the 
stellar lumin osity is proportion al to the halo mass, while the Ly-a luminosity depends on 



Mfjj (x M 2 (IDayal et al.l 120091 ) . The non-linear relation of Ly-a luminosity on the halo 
mass is therefore affected by the scatter in Mi ow /Mhi g h- For example, in our medium-box 
256vs512 simulations at z = 0, (M low /M high ) ~ 0.84 while ((M low /M high ) 2 ) ~ 1.33 for halos 
with 32 < N < 64. This means that the Ly-a luminosity would be overestimated compared 
to the stellar luminosity by a factor ~ 1.6 at low resolution. 



Another example is given by the convergence study carried out in lTrenti et al.l (120091 ) to 
validate their simulations for the transition from metal-free to metal-enriched star formation 
during the reionization epoch. They found that a "low" resolution run (N p = 512 3 ) achieved 
full convergence with the high-resolution run (N p = 1024 3 ) only at z < 7, w hen sub-grid 



physi cs was implemented in DM halos with iV > 100 particles (see Fig. 6 in iTrenti et al. 



20091 ). Individual scatter in halo masses might also introduce numerical noise when feedback 
is considered. For example, if supernova feedback from star formation at the center of the 
halos is near the critical level to evacuate most of the baryons, then the numerical uncertainty 
in the halo mass might play a critical role for the future development of the star formation 
history in the descendant halos of poorly resolved progenitors. 

Overall we recommend implementing extra physics only in halos with iV > 100 for 
rare halos [y > 3) and N > 400 for very common halos [y < 1). This will guarantee 
£(7V") < 0.2 and adequate convergence of other halo properties. If this recommendation is 
followed, then both the Enzo code and Gadget2 are in a regime where their convergence 
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properties are similar. Of course, many interesting problems in cosmology require a greater 
dynamic range than is currently possible to resolve. Thus, implementing extra physics only 
on halos resolved with N > 100 particles is not optimal. In such cases, modelers should take 
precautions to demonstrate convergence of their halo simulations with iV < 100 — 400 and 
avoid extrapolating beyond the range of validity. 

We thank Zarija Lukic, Brian O'Shea and Darren Reed for stimulating comments, Stef- 
fen Knollmann for his advice on compiling and running the Amiga halo finder code, and 
an anonymous referee for useful suggestions. We acknowledge support from the University 
of Colorado Astrophysical Theory Program through grants from NASA (NNX07AG77G) 
and NSF (AST07-07474). This research was supported in part by the National Science 
Foundation through TeraGrid awards TG-AST090040 and TG-AST090101. 
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Table 1: Simulation Summary 
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Note. — Summary of the properties of our cosmological simulations done with Gadget2. The first column 
reports the number of dark matter particles N p , the second the box-size lb ox - The single-particle dark-matter 
mass (judm) is in the third column, the initial redshift z sta rt in the fourth, and the final redshift z en( i in the 
last column. We also carried out the medium-box runs (lbox = 64 Mpc h^ 1 ) with the Enzo code, in both 
Unigrid and AMR (6 levels) mode for runs up to N p = 256 3 and in Unigrid mode only for N p — 512 3 . 
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Fig. 1. — Comparison of individual halo masses at z = for four realizations of the same 
initial conditions (medium box: lf, ox = 64 Mpc h _1 ) at different resolution, from N p = 64 3 to 



N p = 512 3 . The two panels show the ratio Mi ow /M 5 i 2 of halo mass measured at low resolution 
against the value found in the highest resolution run as a function of the halo mass at high 
resolution. Halos have been identified in the high-resolution realization (N p = 512 3 ) and 
matched to counterparts at lower resolution. Left panel shows the median of Mi ow /Mr,i2 as a 
solid line. Dashed red lines delimit the symmetric la region around median. Color coding is 
as follows. Green: N[ ow = 64 3 ; red: Ni ow = 128 3 ; blue: Ni ow = 256 3 . Right panel: zoom into 
the individual M; OU) /M 512 values for the most massive halos (identified in all realizations). 
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Fig. 2. — Upper left panel: Ratio of low-to-high resolution mass of dark matter halos identi- 
fied in the N p = 128 3 version of our medium box ilhox = 64 Mpc h -1 ) at z = and compared 
to higher resolution realization (N p = 512 3 ) with same initial conditions. The scatter plot 
is shown as function of number of low-resolution particles. Solid red line represents median 
value of ratio, and red dashed lines enclose la symmetric confidence region around median. 
The upper axis translates halo particle number into dimensionless variable v = 5 2 /a 2 (M), 
a measure of the rarity of a halo in the extended Press-Schechter formalism. Upper right 
panel: same as upper left panel, but her e dark-matter halos have been identified with Amiga 
halo finder (IKnollmann fc Knebd 120091 ) including removal of unbound particles. Lower left 
and right panels are same as the upper set, but for N p = 256 3 run compared to N p = 512 3 
run. 
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Fig. 3. — Halo mass function at z = of medium box simulated with N p = 512 3 particles 
(black solid line) compared to N p = 256 3 realizations (Gadget2 blue dotted-dashed line; 
Enz o AMR red dotted - dashe d line; Enzo Unigrid cyan dashed line). For comparison, we also 
plot ISheth fc Tormenl ( 119991 ) halo mass function (solid green line). The low- resolution halo 
mass function obtained in the Gadget2 run is consistent with high-resolution and analytical 
mass function down to its resolution limit (N ~ 30). The halo mass function from the Enzo 
runs converges only when N > 400. 
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Fig. 4. — lcr confidence intervals (dashed lines) and median (solid lines) for the ratio of 
low-to-high resolution mass of dark matter halos as in Figure [2] but considering different 
simulations. The left panel investigates different box sizes. Blue lines are associated with 
the N[ ow = 512 3 version of our large box {l hox = 512 Mpc h _1 ) compared to the N high = 
1024 3 realization at z = 0. Red lines refer to Ni ow = 256 3 version of medium box (lf, ox = 
64 Mpc h _1 ) compared to Nhigh = 512 3 realization at z = 0. Green lines refer to Ni ow = 256 3 
version of small box (lb ox = 64 Mpc h _1 ) compared to Nhigh — 512 3 realization at z = 6. 
Right panel investigates the redshift dependence of convergence for large box realization 
Niow = 512 3 compared to Nhigh = 1024 3 at z = (blue), z — 1 (red), and z = 4 (green). 
Halos at fixed number of particles are rarer at higher redshift and have slightly less scatter 
in the mass measurements. 
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Fig. 5. — The la dispersion £ around median for Mi ow /M high distribution in our set of 
simulations considering halos with 100 < Ni ow < 200 (black squares) and with 1000 < 
Ni ow < 2000 (blue pentagons) as function of rarity v. Rare (large v) halos resolved with a 
low number of particles have less scatter than their more common counterparts. 
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Fig. 6. — The la dispersion £ around median for M low /M high distribution in our set of 
simulations as function of number of particles in low-resolution run. Triangles show color- 
coded results from medium box: common halos in the 256vs512 snapshot at z = (green) 
and at z — 1 (black), rarer halos in the z = 4 snapshot (light blue) and halos in the 
128vs512 snapshot at z = (violet). Squares are from large box at z = (red 256vsl024; 
blue 512vsl024). Pentagons are from large box (512vsl024) at z = 1 and blue hexagons 
are from same simulation at z — 4, representing halos with smallest £(iV) (due to their 
rarity). Yellow hexagons are from small box at z = 6 (256vs512). Dashed line represents 
the analytical scaling predicted in Section [3] (£ oc iV -1 / 3 ). 
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Fig. 7. — (7 Mpc /i -1 ) 3 region around one halo of our medium-box simulations at low reso- 
lution (left N p = 256 3 ) and high resolution (right N p = 512 3 ). The 8192 particles belonging 
to the halo at low resolution are shown in magenta (left), while their counterparts at high 
resolution are in red (right). This halo has two distinct counterparts at high resolution, 
because the merging is slightly delayed compared to low-resolution realization, where a thin 
bridge of particles already connects the two sub-components. 
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Fig. 8.— Same as Fig. H] but for runs carried out with the Enzo code. The left panel 
compares Ni ow = 128 3 in Unigrid mode, shown as red lines, and in AMR mode (6 levels), 
shown as blue lines, to Nhigh — 512 3 (Unigrid) at z = 0. The right panel is as the left one but 
for Ni ow = 256 3 . The mass of individual halos with N < 500 is significantly underestimated 
in Enzo. AMR helps to improve convergence, but its effectiveness is limited. 
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Fig. 9. — Convergence of additional properties of dark matter halos in N p = 256 3 realization 
of medium box compared to their high-resolution counterparts (N p = 512 3 run). Upper- 
left panel: Ratio of core density p 2 5 iow/p25 high- Upper-right panel: Ratio of virial radius 
r vir iow/ r vir high- Lower left panel: Ratio of dimensionless angular momentum \i ow /^high- 
For each panel the solid red line shows the median of distribution. Dashed red lines delimit 
the symmetric la region around median. 



